home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / EDUCNOMY / ASTROSET.LZH / COMET.BAS < prev    next >
BASIC Source File  |  1985-12-20  |  6KB  |  236 lines
  1. 5 GOSUB 69: GOTO 200
  2. 10 REM   KEPLER'S EQUATION
  3. 11 REM
  4. 12 P1=3.14159265: R1=180/P1
  5. 13 K=0.01720209895
  6. 18 IF E0>=0.95 THEN GOTO 40
  7. 19 IF E0>=1 THEN 85
  8. 20 A1=Q/(1-E0): M=K*T*A1^(-1.5)
  9. 21 REM
  10. 22 REM   BINARY SEARCH
  11. 23 REM
  12. 24 F=SGN(M): M=ABS(M)/(2*P1)
  13. 25 M=(M-INT(M))*2*P1*F
  14. 26 IF M<0 THEN M=M+2*P1
  15. 27 F=1: IF M>P1 THEN F=-1
  16. 28 IF M>P1 THEN M=2*P1-M
  17. 29 E=P1/2: D=P1/4
  18. 30 FOR I1=1 TO 23
  19. 31 M1=E-E0*SIN(E)
  20. 32 E=E+SGN(M-M1)*D: D=D/2
  21. 33 NEXT I1
  22. 34 V=SQR((1+E0)/(1-E0)): E=E*F
  23. 35 V=2*ATN(V*SIN(E/2)/COS(E/2))
  24. 36 R=A1*(1-E0*COS(E))
  25. 38 GOTO 81
  26. 39 REM
  27. 40 REM    GAUSS METHOD
  28. 41 REM
  29. 43 A=SQR((1+9*E0)/10)
  30. 44 B=5*(1-E0)/(1+9*E0)
  31. 45 C=SQR(5*(1+E0)/(1+9*E0))
  32. 46 B1=3*A*K*T/SQR(2*Q*Q*Q)
  33. 47 B2=1: REM  INITIAL VALUE
  34. 48 W1=B2*B1: B3=ATN(2/W1)
  35. 49 T1=SIN(B3/2)/COS(B3/2)
  36. 50 S1=SGN(T1): T1=ABS(T1)
  37. 51 T2=T1^(1/3)*S1: G=ATN(T2)
  38. 52 S=2*COS(2*G)/SIN(2*G)
  39. 53 A2=B*S*S: B0=B2: B2=0
  40. 54 IF ABS(A2)>0.3 THEN 19
  41. 55 FOR J=0 TO 7
  42. 56 B2=B2+B(J)*A2^J
  43. 57 NEXT J
  44. 58 IF ABS(B2-B0)>1E-8 THEN 48
  45. 59 C1=0
  46. 60 FOR J=0 TO 7
  47. 61 C1=C1+S(J)*A2^J
  48. 62 NEXT J
  49. 63 C1=SQR(1/C1)
  50. 64 V1=C*C1*S: D1=1/(1+A2*C1*C1)
  51. 65 V=2*ATN(V1): R=Q*D1*(1+V1*V1)
  52. 67 GOTO 81
  53. 68 REM
  54. 69 REM   COEFFICIENTS
  55. 70 FOR J=0 TO 7: READ B(J): NEXT
  56. 71 FOR J=0 TO 7: READ S(J): NEXT
  57. 72 RETURN
  58. 73 DATA 1,0,-0.017142857
  59. 74 DATA -0.003809524, -0.001104267
  60. 75 DATA -0.000367358,-0.000131675
  61. 76 DATA -0.000049577,1,-0.8
  62. 77 DATA 0.04571429,0.01523810
  63. 78 DATA 0.00562820, 0.00218783
  64. 79 DATA 0.00087905,0.00036155
  65. 80 RETURN
  66. 81 IF V<0 THEN V=V+2*P1
  67. 84 GOTO 86
  68. 85 PRINT "OUT OF RANGE"
  69. 86 RETURN
  70. 200 REM   COMET POSITION IN SKY
  71. 203 REM
  72. 206 PRINT "PERI DATE  ";
  73. 207 GOSUB 800: J9=J: F9=F
  74. 208 INPUT "PERI DIST Q  ";Q
  75. 209 INPUT "ECCENTRICITY ";E0
  76. 212 INPUT "ARG OF PERIHELION ";W
  77. 215 INPUT "LONG OF ASC NODE  ";N
  78. 218 INPUT "INCLINATION       ";I
  79. 233 P1=3.14159265: R1=P1/180
  80. 236 E=23.4457889*R1: REM OBLIQUITY
  81. 239 W=W*R1: N=N*R1: I=I*R1
  82. 242 GOSUB 338  
  83. 245 PRINT
  84. 246 PRINT "DATE  ";
  85. 247 GOSUB 800: J1=J: F1=F
  86. 248 GOSUB 500: T=(J1-J9)+(F1-F9)
  87. 249 GOSUB 10
  88. 251 REM   POSITION IN ORBIT PLANE
  89. 254 X1=R*COS(V): Y1=R*SIN(V)
  90. 257 REM
  91. 260 REM   HELIOCENTRIC EQUATORIAL
  92. 261 REM   COORDINATES
  93. 263 X2=P7*X1+Q7*Y1
  94. 266 Y2=P8*X1+Q8*Y1
  95. 269 Z2=P9*X1+Q9*Y1
  96. 272 REM
  97. 275 REM   GEOCENTRIC EQUATORIAL
  98. 276 REM   COORDINATES
  99. 278 X3=X+X2: Y3=Y+Y2: Z3=Z+Z2
  100. 281 GOSUB 392: REM  FOR 2000.0!
  101. 284 D3=SQR(X3*X3+Y3*Y3+Z3*Z3)
  102. 287 R1=P1/180
  103. 290 A=ATN(Y3/X3)/(15*R1)
  104. 293 IF X3<0 THEN A=A+12
  105. 296 IF A<0 THEN A=A+24
  106. 299 D=ATN(Z3/SQR(X3*X3+Y3*Y3))/R1
  107. 301 REM   NOW ROUND OFF
  108. 302 A=A+0.05/60
  109. 305 H=INT(A): M=60*(A-H)
  110. 308 M=INT(10*M)/10
  111. 311 S=SGN(D): D=ABS(D)+0.5/60
  112. 314 D1=INT(D): M1=INT(60*(D-D1))
  113. 317 S$="+": IF S=-1 THEN S$="-"
  114. 320 PRINT "R.A.  ";H;"H ";M;"M"
  115. 323 PRINT "DEC. ";S$;D1;"D ";M1;"M"
  116. 326 PRINT "COMET-EARTH DIST:";D3
  117. 329 PRINT "COMET-SUN DIST: ";R
  118. 330 INPUT "ANOTHER (Y OR N) ";Q$
  119. 331 IF Q$="Y" THEN 245
  120. 332 END
  121. 335 REM
  122. 338 REM   P'S AND Q'S
  123. 344 W1=SIN(W): W2=COS(W)
  124. 347 N1=SIN(N): N2=COS(N)
  125. 350 I1=SIN(I): I2=COS(I)
  126. 353 E1=SIN(E): E2=COS(E)
  127. 356 P7=W2*N2-W1*N1*I2
  128. 359 P8=(W2*N1+W1*N2*I2)*E2
  129. 362 P8=P8-W1*I1*E1
  130. 365 P9=(W2*N1+W1*N2*I2)*E1
  131. 368 P9=P9+W1*I1*E2
  132. 371 Q7=-W1*N2-W2*N1*I2
  133. 374 Q8=(-W1*N1+W2*N2*I2)*E2
  134. 377 Q8=Q8-W2*I1*E1
  135. 380 Q9=(-W1*N1+W2*N2*I2)*E1
  136. 383 Q9=Q9+W2*I1*E2
  137. 386 RETURN
  138. 389 REM
  139. 392 REM   1950.0 --> 2000.0
  140. 395 A7=+0.9999257: A8=-0.0111789
  141. 398 A9=-0.0048590: B7=+0.0111789
  142. 401 B8=+0.9999375: B9=-0.0000272
  143. 404 C7=+0.0048590: C8=-0.0000272
  144. 407 C9=+0.9999882
  145. 410 X4=A7*X3+A8*Y3+A9*Z3
  146. 413 Y4=B7*X3+B8*Y3+B9*Z3
  147. 416 Z4=C7*X3+C8*Y3+C9*Z3
  148. 419 X3=X4: Y3=Y4: Z3=Z4
  149. 422 RETURN
  150. 500 REM    X,Y,Z OF THE SUN
  151. 501 REM    (EQUINOX 1950.0)
  152. 502 REM
  153. 504 J8=J-2415020: R1=3.14159265/180
  154. 505 T=(J8+F)/36525
  155. 506 P0=1.396041+0.000308*(T+0.5)
  156. 507 P0=P0*(T-0.499998)
  157. 508 A=100:GOSUB 529:G0=A+358.475833
  158. 509 L0=A+279.696678-P0
  159. 510 A=1336: GOSUB 529  
  160. 511 C0=A+270.434164-P0
  161. 512 A=162: GOSUB 529  
  162. 513 V0=A+212.603219
  163. 514 A=53:GOSUB 529: M0=A+319.529425
  164. 515 A=8: GOSUB 529: J0=A+225.444651
  165. 516 G=G0+T*(-0.950250-0.000150*T)
  166. 517 C=C0+T*(307.883142-0.001133*T)
  167. 518 L=L0+T*(0.768920+0.000303*T)
  168. 519 V=V0+T*(197.803875+0.001286*T)
  169. 520 M=M0+T*(59.8585+0.000181*T)
  170. 521 J=J0+T*154.906654
  171. 522 G=G*R1: C=C*R1: L=L*R1
  172. 523 V=V*R1: M=M*R1: J=J*R1
  173. 524 GOSUB 532  
  174. 528 RETURN
  175. 529 REM   NORMALIZATION
  176. 530 A=360*(A*T-INT(A*T)): RETURN
  177. 531 REM
  178. 532 X=0.000011*COS(2*G-L-2*J)
  179. 533 X=X+0.000011*COS(2*G+L-2*V)
  180. 534 X=X-0.000012*COS(G+L-V)
  181. 535 X=X-0.000012*COS(4*G-L-8*M+3*J)
  182. 536 X=X+0.000012*COS(4*G+L-8*M+3*J)
  183. 537 X=X-0.000014*COS(C-2*L)
  184. 538 X=X+0.000017*COS(C)
  185. 539 X=X+0.000018*SIN(2*G+L-2*V)
  186. 540 X=X-0.000021*T*COS(G+L)
  187. 541 X=X-0.000026*SIN(G-L-J)
  188. 542 X=X+0.000035*COS(2*G-L)
  189. 543 X=X+0.000063*T*COS(G-L)
  190. 544 X=X+0.000105*COS(2*G+L)
  191. 545 X=X+0.008374*COS(G+L)
  192. 546 X=X-0.025127*COS(G-L)
  193. 547 X=X+0.999860*COS(L)
  194. 548 REM
  195. 549 Y=0.000010*SIN(2*G+L-2*V)
  196. 550 Y=Y-0.000010*SIN(2*G-L-2*J)
  197. 551 Y=Y-0.000011*SIN(G+L-V)
  198. 552 Y=Y+0.000011*SIN(4*G-L-8*M+3*J)
  199. 553 Y=Y+0.000011*SIN(4*G+L-8*M+3*J)
  200. 554 Y=Y+0.000013*SIN(C-2*L)
  201. 555 Y=Y+0.000016*SIN(C)
  202. 556 Y=Y-0.000017*COS(2*G+L-2*V)
  203. 557 Y=Y-0.000019*T*SIN(G+L)
  204. 558 Y=Y-0.000024*COS(G-L-J)
  205. 559 Y=Y-0.000032*SIN(2*G-L)
  206. 560 Y=Y-0.000057*T*SIN(G-L)
  207. 561 Y=Y+0.000097*SIN(2*G+L)
  208. 562 Y=Y+0.007683*SIN(G+L)
  209. 563 Y=Y+0.023053*SIN(G-L)
  210. 564 Y=Y+0.917308*SIN(L)
  211. 565 REM
  212. 566 Z=-0.000010*COS(G-L-J)
  213. 567 Z=Z-0.000014*SIN(2*G-L)
  214. 568 Z=Z-0.000025*T*SIN(G-L)
  215. 569 Z=Z+0.000042*SIN(2*G+L)
  216. 570 Z=Z+0.003332*SIN(G+L)
  217. 571 Z=Z+0.009998*SIN(G-L)
  218. 572 Z=Z+0.397825*SIN(L)
  219. 573 RETURN
  220. 800 REM   CALENDAR --> JD
  221. 805 REM
  222. 810 INPUT "Y,M,D ";Y,M,D
  223. 815 G=1
  224. 820 D1=INT(D): F=D-D1-0.5
  225. 825 J=-INT(7*(INT((M+9)/12)+Y)/4)
  226. 830 IF G=0 THEN 850  
  227. 835 S=SGN(M-9): A=ABS(M-9)
  228. 840 J3=INT(Y+S*INT(A/7))
  229. 845 J3=-INT((INT(J3/100)+1)*3/4)
  230. 850 J=J+INT(275*M/9)+D1+G*J3
  231. 855 J=J+1721027+2*G+367*Y
  232. 860 IF F>=0 THEN 870  
  233. 865 F=F+1: J=J-1
  234. 870 RETURN
  235. 875 END
  236.